Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation.
Undertake the analysis of the primary outcome for ASCOT.
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
theme_set(theme_minimal(base_size = 12))
bayesplot_theme_set(theme_minimal(base_size = 12))
color_scheme_set("red")
options(digits = 4)source("r/derive_full_datasets.r")
all_dat <- read_all_no_daily()The primary outcome is a composite of death, need for new respiratory support, or vasopressor/inotropic support. From the protocol:
Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation. This includes any participant who receives non-invasive mechanical ventilation (either CPAP or BIPAP, apart from the below considerations) any time after enrolment even if not transferred to ICU. It does NOT include the use of humidified high-flow nasal prong oxygen.
Participants on pre-existing home BiPAP or CPAP will not be considered to have met the primary outcome unless they have either:
- required invasive mechanical ventilation (i.e. intubation), or
- graduated from CPAP only whilst asleep to BiPAP at any time, or
- graduated from BiPAP only whilst asleep to BiPAP for >12 hours/day, or
- died by day 28
There may be cases where a patient has been assessed as requiring intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support, but the patient or family declined treatment and the patient was discharged home. If attempts to obtain 28-day data are unsuccessful or not possible, and the investigator had deemed at the time of discharge that the patient would be highly likely to die within 28 days from randomisation, these participants will be deemed to have met the primary outcome.
Derivation of the outcome requires checking of the daily, discharge, and day 28 data extracts. On the daily data, there is a variable DD_PrimaryEndpointReachedToday however this was coded incorrectly in the original database and therefore fails to capture some participants. Additionally, given the composite nature of the outcome it is useful to check all components individually as well as the composite outcome. Therefore, this variable is not used in the derivation, but is included as a cross-check.
Below, each component is summarised in aggregate.
There were 50 participants who had died by day 28, 40 were are reported as a discharge status of death and 10 were reported as a death post-discharge. Day 28 mortality was missing or unknown for 30 participants.
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
select(DIS_death, DIS_day, D28_death) %>%
dplyr::count(DIS_death, DIS_day <= 28, D28_death) %>%
spread(D28_death, n, fill = 0)# A tibble: 4 × 5
DIS_death `DIS_day <= 28` `0` `1` `<NA>`
<int> <lgl> <dbl> <dbl> <dbl>
1 0 FALSE 10 0 0
2 0 TRUE 1489 10 30
3 1 FALSE 3 0 0
4 1 TRUE 0 40 0Vasopressor/inotropic support is reported on both the daily (from day 1 to discharge) and the day 28 (from discharge to day 28) extracts. There were 39 participants who required vasopressor/inotropic support. Of these, 19 were reported between discharge and day 28, 18 were reported prior to discharge, and 2 were reported for both. The component was missing for 35 participants. Either due to missing information at day 28 (33) or due to missing daily information (2).
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
select(DAILY_missing, ANY_vasop, ANY_DD_vasop, D28_vasop) %>%
dplyr::count(DAILY_missing, ANY_DD_vasop, D28_vasop, ANY_vasop) %>%
spread(ANY_vasop, n, fill = 0)# A tibble: 7 × 6
DAILY_missing ANY_DD_vasop D28_vasop `0` `1` `<NA>`
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 0 1508 0 0
2 0 0 1 0 19 0
3 0 0 NA 0 0 33
4 0 1 0 0 18 0
5 0 1 1 0 2 0
6 1 NA 0 0 0 1
7 1 NA NA 0 0 1Respiratory support is reported on both the daily (from day 1 to discharge) and the day 28 (at day 28) extracts. There were 70 participants who required new intensive respiratory support. Of these, 60 were prior to discharge, 6 were reported at day 28, and 4 were both. The outcome was unknown for 36 participants due to either missing day 28 information (34) or missing daily information (2).
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
select(DAILY_missing, ANY_vent, ANY_DD_vent, D28_vent) %>%
dplyr::count(DAILY_missing, ANY_DD_vent, D28_vent, ANY_vent) %>%
spread(ANY_vent, n, fill = 0)# A tibble: 6 × 6
DAILY_missing ANY_DD_vent D28_vent `0` `1` `<NA>`
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 0 1476 0 0
2 0 0 NA 0 0 34
3 0 1 0 0 60 0
4 0 1 1 0 4 0
5 0 1 NA 0 6 0
6 1 NA 0 0 0 2There were 3 participants who were DAMA and identified as likely to die. However, the day 28 status was observed for all 3 participants, therefore, no-one met the primary endpoint due to DAMA and unknown day 28 status.
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
select(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
dplyr::count(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
spread(D28_death, n)# A tibble: 3 × 5
DIS_DAMA DIS_DAMAlikelytodie `0` `1` `<NA>`
<int> <int> <int> <int> <int>
1 0 0 1469 46 28
2 1 0 32 2 2
3 1 1 1 2 NAThere were 107 (7% of observed outcomes) participants who met the composite outcome and 36 (2% of all outcomes) participants whose outcome was unknown. Of the 107 who met the outcome, 12 met all three criteria (death, vasopressor, ventilation) and the remaining met some subset of them. Of the 36 with missing outcome, 26 were missing all three components.
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(PO)# A tibble: 3 × 2
PO n
<dbl> <int>
1 0 1439
2 1 107
3 NA 36all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
select(PO, D28_death, ANY_vasop, ANY_vent) %>%
dplyr::count(D28_death, ANY_vasop, ANY_vent, PO) %>%
spread(PO, n, fill = 0)# A tibble: 15 × 6
D28_death ANY_vasop ANY_vent `0` `1` `<NA>`
<int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 0 1439 0 0
2 0 0 1 0 31 0
3 0 1 0 0 20 0
4 0 1 1 0 5 0
5 0 NA 0 0 0 5
6 0 NA NA 0 0 2
7 1 0 0 0 10 0
8 1 0 1 0 21 0
9 1 0 NA 0 4 0
10 1 1 0 0 2 0
11 1 1 1 0 12 0
12 1 NA NA 0 1 0
13 NA 0 NA 0 0 3
14 NA NA 1 0 1 0
15 NA NA NA 0 0 26agedat <- all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(PO, AgeAtEntry) %>%
spread(PO, n, fill = 0) %>%
mutate(p = `1` / (`1` + `0`))
agemod <- glm(
cbind(`1`, `0`) ~ AgeAtEntry,
data = agedat,
family = binomial())
agedat <- agedat %>%
mutate(
ypred = predict(agemod, newdata = agedat, type = "response")
)
p <- ggplot(agedat, aes(AgeAtEntry, p)) +
geom_point() +
geom_line(aes(y = ypred)) +
labs(y = "Proportion with outcome", x = "Age at entry")
ggplotly(p)all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(`Age 60+` = agegte60, PO) %>%
group_by(`Age 60+`) %>%
spread(PO, n, fill = 0) %>%
mutate(
p_1 = `1` / (`1` + `0`),
p_miss = `<NA>` / (`1` + `0` + `<NA>`)
) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped")| Age 60+ | 0 | 1 | <NA> | p_1 | p_miss |
|---|---|---|---|---|---|
| 0 | 1038 | 59 | 25 | 0.05 | 0.02 |
| 1 | 401 | 48 | 11 | 0.11 | 0.02 |
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(Country = PT_CountryName, PO) %>%
group_by(Country) %>%
spread(PO, n, fill = 0) %>%
mutate(
p_1 = `1` / (`1` + `0`),
p_miss = `<NA>` / (`1` + `0` + `<NA>`)
) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped")| Country | 0 | 1 | <NA> | p_1 | p_miss |
|---|---|---|---|---|---|
| Australia | 128 | 8 | 14 | 0.06 | 0.09 |
| India | 1190 | 83 | 15 | 0.07 | 0.01 |
| Nepal | 104 | 11 | 3 | 0.10 | 0.03 |
| New Zealand | 17 | 5 | 4 | 0.23 | 0.15 |
tabdat <- all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(
Country = PT_CountryName,
Site = PT_LocationName,
PO) %>%
group_by(Country, Site) %>%
spread(PO, n, fill = 0) %>%
mutate(
p_1 = `1` / (`1` + `0`),
p_miss = `<NA>` / (`1` + `0` + `<NA>`)
) %>%
ungroup()
n_country <- tabdat %>%
dplyr::count(Country) %>%
pull(n)
cn <- cumsum(n_country)
tabdat %>%
select(-Country) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped") %>%
pack_rows("Australia", 1, cn[1]) %>%
pack_rows("India", cn[1] + 1, cn[2]) %>%
pack_rows("Nepal", cn[2] + 1, cn[3]) %>%
pack_rows("New Zealand", cn[3] + 1, cn[4])| Site | 0 | 1 | <NA> | p_1 | p_miss |
|---|---|---|---|---|---|
| Australia | |||||
| Alfred Hospital | 6 | 1 | 7 | 0.14 | 0.50 |
| Blacktown Hospital | 4 | 0 | 0 | 0.00 | 0.00 |
| Box Hill Hospital | 4 | 2 | 2 | 0.33 | 0.25 |
| Campbelltown Hospital | 11 | 1 | 0 | 0.08 | 0.00 |
| John Hunter Hospital | 2 | 0 | 0 | 0.00 | 0.00 |
| Liverpool Hospital | 6 | 2 | 0 | 0.25 | 0.00 |
| Monash Health | 8 | 0 | 0 | 0.00 | 0.00 |
| Prince Charles Hospital | 3 | 0 | 0 | 0.00 | 0.00 |
| Royal Melbourne Hospital | 7 | 1 | 1 | 0.12 | 0.11 |
| Royal North Shore Hospital | 6 | 0 | 0 | 0.00 | 0.00 |
| Royal Perth Hospital | 5 | 0 | 1 | 0.00 | 0.17 |
| St George Hospital | 1 | 0 | 0 | 0.00 | 0.00 |
| St Vincent's Hospital, Sydney | 10 | 0 | 1 | 0.00 | 0.09 |
| Wagga Wagga Base Hospital | 10 | 0 | 0 | 0.00 | 0.00 |
| Western Health | 1 | 0 | 0 | 0.00 | 0.00 |
| Westmead Hospital | 42 | 1 | 2 | 0.02 | 0.04 |
| Wollongong Hospital | 2 | 0 | 0 | 0.00 | 0.00 |
| India | |||||
| Aditya Multispeciality Hospital | 63 | 0 | 0 | 0.00 | 0.00 |
| Apex Hospitals, Malviya Nagar | 125 | 1 | 0 | 0.01 | 0.00 |
| Believers Church Medical College Hospital | 18 | 0 | 0 | 0.00 | 0.00 |
| Christian Medical College Vellore | 84 | 16 | 10 | 0.16 | 0.09 |
| CMC Ludhiana | 64 | 18 | 0 | 0.22 | 0.00 |
| Core Hospital | 86 | 21 | 0 | 0.20 | 0.00 |
| Jivanrekha Multispeciality Hospital | 329 | 13 | 3 | 0.04 | 0.01 |
| Maharaja Agrasen Superspeciality Hospital | 138 | 0 | 0 | 0.00 | 0.00 |
| Samishta Hospital & Research Institute | 156 | 2 | 0 | 0.01 | 0.00 |
| Sterling Multispecialty Hospital Pune | 85 | 12 | 2 | 0.12 | 0.02 |
| Symbiosis University Hospital | 42 | 0 | 0 | 0.00 | 0.00 |
| Nepal | |||||
| Bir Hospital | 71 | 8 | 3 | 0.10 | 0.04 |
| Tribhuvan University Teaching Hospital | 33 | 3 | 0 | 0.08 | 0.00 |
| New Zealand | |||||
| Christchurch Hospital | 1 | 0 | 0 | 0.00 | 0.00 |
| Dunedin Hospital | 1 | 0 | 1 | 0.00 | 0.50 |
| Middlemore Hospital | 13 | 5 | 0 | 0.28 | 0.00 |
| North Shore Hospital | 2 | 0 | 3 | 0.00 | 0.60 |
caldat <- all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(PO, yr = year(RandDate), mth = month(RandDate)) %>%
spread(PO, n, fill = 0) %>%
mutate(p = `1` / (`1` + `0`))
p <- ggplot(caldat, aes(mth, p)) +
facet_wrap( ~ yr) +
geom_point() +
labs(
y = "Proportion with outcome",
x = "Calendar date (year and month)") +
scale_x_continuous(breaks = 1:12)
ggplotly(p)all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(`Anti-coagulation` = CAssignment, PO) %>%
group_by(`Anti-coagulation`) %>%
spread(PO, n, fill = 0) %>%
mutate(
p_1 = `1` / (`1` + `0`),
p_miss = `<NA>` / (`1` + `0` + `<NA>`)
) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped")| Anti-coagulation | 0 | 1 | <NA> | p_1 | p_miss |
|---|---|---|---|---|---|
| C0 | 16 | 4 | 6 | 0.20 | 0.23 |
| C1 | 551 | 45 | 14 | 0.08 | 0.02 |
| C2 | 571 | 30 | 12 | 0.05 | 0.02 |
| C3 | 258 | 21 | 4 | 0.08 | 0.01 |
| C4 | 43 | 7 | 0 | 0.14 | 0.00 |
all_dat %>%
filter(ENR_rec == 1, WTH_FU == 0) %>%
dplyr::count(
Country = PT_CountryName,
`Anti-coagulation` = CAssignment,
PO) %>%
complete(Country, `Anti-coagulation`, PO, fill = list(n = 0)) %>%
group_by(Country, `Anti-coagulation`) %>%
spread(PO, n, fill = 0) %>%
ungroup() %>%
mutate(
p_1 = `1` / (`1` + `0`),
p_miss = `<NA>` / (`1` + `0` + `<NA>`)
) %>%
select(-Country) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped") %>%
pack_rows("Australia", 1, 5) %>%
pack_rows("India", 6, 10) %>%
pack_rows("Nepal", 11, 15) %>%
pack_rows("New Zealand", 16, 20)| Anti-coagulation | 0 | 1 | <NA> | p_1 | p_miss |
|---|---|---|---|---|---|
| Australia | |||||
| C0 | 16 | 2 | 6 | 0.11 | 0.25 |
| C1 | 45 | 2 | 2 | 0.04 | 0.04 |
| C2 | 50 | 3 | 6 | 0.06 | 0.10 |
| C3 | 7 | 0 | 0 | 0.00 | 0.00 |
| C4 | 10 | 1 | 0 | 0.09 | 0.00 |
| India | |||||
| C0 | 0 | 0 | 0 | NaN | NaN |
| C1 | 449 | 37 | 7 | 0.08 | 0.01 |
| C2 | 486 | 25 | 5 | 0.05 | 0.01 |
| C3 | 251 | 21 | 3 | 0.08 | 0.01 |
| C4 | 4 | 0 | 0 | 0.00 | 0.00 |
| Nepal | |||||
| C0 | 0 | 0 | 0 | NaN | NaN |
| C1 | 49 | 4 | 3 | 0.08 | 0.05 |
| C2 | 30 | 1 | 0 | 0.03 | 0.00 |
| C3 | 0 | 0 | 0 | NaN | NaN |
| C4 | 25 | 6 | 0 | 0.19 | 0.00 |
| New Zealand | |||||
| C0 | 0 | 2 | 0 | 1.00 | 0.00 |
| C1 | 8 | 2 | 2 | 0.20 | 0.17 |
| C2 | 5 | 1 | 1 | 0.17 | 0.14 |
| C3 | 0 | 0 | 1 | NaN | 1.00 |
| C4 | 4 | 0 | 0 | 0.00 | 0.00 |
In what follows, the model is built up from the simplest (treatment-only) to the most complex (primary model).
seed <- 17250
logistic_mod <- cmdstan_model("stan/logistic_bernoulli.stan")
make_model_data_trt_only <- function(dat) {
X <- model.matrix(
~ CAssignment,
data = dat[!is.na(dat$PO), ],
contrast = list(CAssignment = contr.orthonorm)
)
y <- as.integer(dat$PO)
y <- y[!is.na(y)]
N <- dim(X)[1]
K <- dim(X)[2]
list(N = N, K = K, X = X, y = y, beta_sd = c(2.5, 1, 1, 1))
}
mod_dat <- make_model_data_trt_only(
all_dat %>%
filter(ENR_rec == 1, ACS_ITT == 1, WTH_FU == 0)
)
snk <- capture.output(
trt_only_fit <- logistic_mod$sample(
data = mod_dat,
seed = seed,
refresh = 0,
iter_sampling = 2500,
chains = 8,
parallel_chains = 8
))Diagnostics
num_divergent num_max_treedepth ebfmi
[1,] 0 0 1.073
[2,] 0 0 1.088
[3,] 0 0 1.029
[4,] 0 0 1.125
[5,] 0 0 1.000
[6,] 0 0 1.054
[7,] 0 0 1.103
[8,] 0 0 1.096trt_only_fit$summary("beta")# A tibble: 4 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 beta[1] -2.48 -2.48 0.132 0.131 -2.70 -2.27 1.00 12751. 13864.
2 beta[2] 0.574 0.582 0.344 0.343 -0.00605 1.12 1.00 13362. 13385.
3 beta[3] -0.0352 -0.0360 0.196 0.198 -0.357 0.287 1.00 13470. 14218.
4 beta[4] -0.413 -0.413 0.209 0.209 -0.756 -0.0712 1.00 14124. 13591.mcmc_trace(trt_only_fit$draws("beta_raw", format = "matrix"))
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(trt_only_fit$draws("beta"))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
transform_coding(C, cbind(1, contr.treatment(4))) %**%
modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))mcmc_hist(modrv["OR"], binwidth = 0.02) +
labs(x = "Odds ratio (vs C1)") +
geom_vline(xintercept = 1, linetype = 2) +
scale_x_log10()
Under this model a component is added for country (region), \(\mathbb{I}[\text{age}\geq 60]\) (sample-mean centred), and intervention specific ineligibility (C3 only)
New terms:
\[ \begin{aligned} \rho_1 &= 0 \\ \rho_r &\overset{\text{iid}}{\sim} \text{Normal}(0, 1), r=2,...,R,\\ \beta_{\text{age}\geq 60} &\sim \text{Normal}(0, 2.5^2) \\ \beta_{\text{inelgC3}} &\sim \text{Normal}(0, 10) \end{aligned} \]
seed <- 71236
logistic_mod <- cmdstan_model("stan/logistic_bernoulli.stan")
make_model_data_noepochsite <- function(dat) {
X <- model.matrix(
~ CAssignment + agegte60_c + country + inelgc3,
data = dat[!is.na(dat$PO), ],
contrast = list(CAssignment = contr.orthonorm)
)
y <- as.integer(dat$PO[!is.na(dat$PO)])
N <- dim(X)[1]
K <- dim(X)[2]
list(N = N, K = K, X = X, y = y,
beta_sd = c(2.5, 1, 1, 1, 2.5, 1, 1, 1, 10))
}
mod_dat <- make_model_data_noepochsite(
all_dat %>%
filter(ENR_rec == 1, ACS_ITT == 1, WTH_FU == 0) %>%
mutate(
age_c = AgeAtEntry - mean(AgeAtEntry),
agegte60_c = agegte60 - mean(agegte60),
country = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
inelgc3 = if_else(EL_inelg_c3 == 1 | is.na(EL_inelg_c3), 1, 0)
)
)
snk <- capture.output(
no_epoch_site_fit <- logistic_mod$sample(
data = mod_dat,
seed = seed,
refresh = 0,
iter_sampling = 2500,
chains = 8,
parallel_chains = 8
))no_epoch_site_fit$summary("beta")# A tibble: 9 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 beta[1] -2.49 -2.48 0.181 0.180 -2.79 -2.19 1.00 12045. 13850.
2 beta[2] 0.547 0.550 0.403 0.402 -0.122 1.20 1.00 13742. 13325.
3 beta[3] 0.0341 0.0338 0.202 0.202 -0.294 0.370 1.00 16197. 15653.
4 beta[4] -0.323 -0.324 0.212 0.211 -0.669 0.0217 1.00 21503. 16268.
5 beta[5] 0.760 0.759 0.209 0.209 0.417 1.10 1.00 22417. 15132.
6 beta[6] -0.282 -0.269 0.407 0.408 -0.974 0.358 1.00 21069. 14367.
7 beta[7] 0.248 0.253 0.391 0.388 -0.402 0.883 1.00 16499. 15453.
8 beta[8] 0.573 0.601 0.597 0.587 -0.447 1.51 1.00 21254. 13934.
9 beta[9] -0.304 -0.297 0.260 0.258 -0.742 0.116 1.00 16008. 15128.Diagnostics
num_divergent num_max_treedepth ebfmi
[1,] 0 0 1.0572
[2,] 0 0 1.0200
[3,] 0 0 1.0122
[4,] 0 0 0.9957
[5,] 0 0 1.1094
[6,] 0 0 1.0026
[7,] 0 0 1.0076
[8,] 0 0 1.0577mcmc_trace(no_epoch_site_fit$draws("beta_raw", format = "matrix"))
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(no_epoch_site_fit$draws("beta"))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
transform_coding(C, cbind(1, contr.treatment(4))) %**%
modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))mcmc_hist(modrv["OR"], binwidth = 0.02) +
labs(x = "Odds ratio (vs C1)") +
geom_vline(xintercept = 1, linetype = 2) +
scale_x_log10()
tibble(
Parameter = c("OR[C2]", "OR[C3]", "OR[C4]"),
`Pr(OR < 1)` = Pr(modrv$OR < 1),
Median = median(modrv$OR),
Mean = E(modrv$OR),
`95% CrI` = apply(
quantile(modrv$OR, c(0.025, 0.975)), 2,
function(x) sprintf("(%.2f, %.2f)", x[1], x[2]))
) %>%
kable(digits = 2, align = "lrrrr") %>%
kable_styling(bootstrap_options = "striped")| Parameter | Pr(OR < 1) | Median | Mean | 95% CrI |
|---|---|---|---|---|
| OR[C2] | 0.96 | 0.65 | 0.67 | (0.40, 1.04) |
| OR[C3] | 0.68 | 0.88 | 0.91 | (0.50, 1.53) |
| OR[C4] | 0.17 | 1.58 | 1.75 | (0.60, 3.92) |
Under this model a component is added for site nested within country.
New terms
\[ \begin{aligned} \xi_{rs} &\overset{\text{iid}}{\sim}\text{Normal}(0,\sigma^2_{r}), s=1,...,S_r, r=1,...,R \\ \sigma_r &\overset{\text{iid}}{\sim} \text{Half-}t(3, 1). \end{aligned} \]
seed <- 71236
noepoch_mod <- cmdstan_model("stan/primary_logistic_model_site.stan")
make_model_data_noepoch <- function(dat) {
site_counts <- dat %>%
dplyr::count(
Country = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
Location)
merge_aus <- site_counts %>%
filter(Country == "AU", n < 5) %>%
pull(Location)
merge_nz <- site_counts %>%
filter(Country == "NZ", n < 5) %>%
pull(Location)
dat <- dat %>%
mutate(
ctry = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
# group sites with < 5 counts
site = fct_collapse(
factor(Location, levels = site_counts$Location),
`AU other` = merge_aus,
`NZ other` = merge_nz
)
)
region_site <- dat %>%
dplyr::count(ctry, site) %>%
mutate(
ctry_num = as.numeric(ctry),
site_num = as.numeric(fct_inorder(site))
)
dat <- dat %>%
left_join(region_site %>% select(-n), by = c("ctry", "site"))
region <- dat$ctry_num
M_region <- max(region)
site <- dat$site_num
M_site <- max(site)
region_by_site <- region_site$ctry_num
X <- model.matrix(
~ CAssignment + agegte60_c + inelgc3,
data = dat,
contrast = list(CAssignment = contr.orthonorm)
)
y <- as.integer(dat$PO)
N <- dim(X)[1]
K <- dim(X)[2]
list(N = N, K = K, X = X, y = y,
M_region = M_region, region = region,
M_site = M_site, site = site,
region_by_site = region_by_site,
beta_sd = c(2.5, 1, 1, 1, 2.5, 10))
}
sub_dat <- all_dat %>%
filter(
ENR_rec == 1,
ACS_ITT == 1,
WTH_FU == 0,
!is.na(PO)) %>%
mutate(
age_c = AgeAtEntry - mean(AgeAtEntry),
agegte60_c = agegte60 - mean(agegte60),
country = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
inelgc3 = if_else(EL_inelg_c3 == 1 | is.na(EL_inelg_c3), 1, 0)
)
mod_dat <- make_model_data_noepoch(sub_dat)
snk <- capture.output(
no_epoch_fit <- noepoch_mod$sample(
data = mod_dat,
seed = seed,
refresh = 0,
iter_sampling = 2500,
chains = 8,
parallel_chains = 8,
adapt_delta = 0.9
))Diagnostics
num_divergent num_max_treedepth ebfmi
[1,] 0 0 0.8572
[2,] 0 0 0.7903
[3,] 0 0 0.8450
[4,] 0 0 0.7790
[5,] 0 0 0.8532
[6,] 0 0 0.7963
[7,] 0 0 0.7857
[8,] 0 0 0.8210no_epoch_fit$summary("beta")# A tibble: 6 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 beta[1] -3.37 -3.36 0.521 0.515 -4.23 -2.54 1.00 5013. 7718.
2 beta[2] 0.680 0.684 0.428 0.430 -0.0357 1.38 1.00 19585. 15465.
3 beta[3] 0.106 0.105 0.214 0.216 -0.245 0.457 1.00 18640. 15799.
4 beta[4] -0.258 -0.256 0.228 0.228 -0.631 0.116 1.00 22127. 14602.
5 beta[5] 0.682 0.684 0.229 0.231 0.306 1.06 1.00 25053. 14777.
6 beta[6] 0.478 0.483 0.309 0.308 -0.0355 0.978 1.00 18180. 15850.mcmc_trace(no_epoch_fit$draws(
c("beta_raw", "beta_region_raw"),
format = "matrix"))
mcmc_trace(no_epoch_fit$draws(
c("epsilon_site", "tau_site"),
format = "matrix"))
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(
no_epoch_fit$draws(c("beta", "beta_region", "gamma_site", "tau_site")))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
transform_coding(C, cbind(1, contr.treatment(4))) %**%
modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))mcmc_hist(modrv["OR"], binwidth = 0.02, freq = F) +
labs(x = "Odds ratio (vs C1)") +
geom_vline(xintercept = 1, linetype = 2) +
scale_x_log10()
tibble(
Parameter = c("OR[C2]", "OR[C3]", "OR[C4]"),
`Pr(OR < 1)` = Pr(modrv$OR < 1),
Median = median(modrv$OR),
Mean = E(modrv$OR),
`95% CrI` = apply(
quantile(modrv$OR, c(0.025, 0.975)), 2,
function(x) sprintf("(%.2f, %.2f)", x[1], x[2]))
) %>%
kable(digits = 2, align = "lrrrr") %>%
kable_styling(bootstrap_options = "striped")| Parameter | Pr(OR < 1) | Median | Mean | 95% CrI |
|---|---|---|---|---|
| OR[C2] | 0.98 | 0.61 | 0.62 | (0.37, 0.99) |
| OR[C3] | 0.85 | 0.72 | 0.76 | (0.40, 1.30) |
| OR[C4] | 0.18 | 1.60 | 1.81 | (0.57, 4.21) |